Evaluating temperature-dependence of growth allometry
Max Lindmark, Jan Ohlberger, Anna GĂ„rdmark
15 November, 2021
Intro
Fit log-linear models to growth in line with MTE but with random effects, mass-temperature interactions and within species data.
# Load libraries, install if needed
library(tidyverse)
library(brms)
library(RCurl)
library(tidybayes)
library(bayesplot)
library(RColorBrewer)
library(viridis)
library(tidylog)
library(modelr)
library(patchwork)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
options(mc.cores = parallel::detectCores())
# To load entire cache in interactive r session, do: qwraps2::lazyload_cache_dir(path = "R/allometric_growth_model/html")Read and explore data
- decide what to do with data beyond peak for a size groupâŠ
# Read in your data file
dat <-
read.csv(text = getURL("https://raw.githubusercontent.com/maxlindmark/scaling/master/data/growth_analysis.csv"))
dat <- dat %>%
group_by(species_ab) %>%
filter(y > 0) %>%
mutate(log_y = log(y),
log_mass_intra = log_mass - mean(log_mass),
temp_arr_intra = temp_arr - mean(temp_arr)) %>%
ungroup()
dat <- dat %>%
filter(y > 0) %>%
mutate(log_y = log(y),
log_mass_ct = log_mass - mean(log_mass),
temp_arr_ct = temp_arr - mean(temp_arr)) %>%
ungroup()
# Filter data points at below optimum temperatures
#dat <- dat %>% filter(above_peak_temp == "N")
colnames(dat)
#> [1] "y" "growth_rate_..day" "geom_mean_mass_g"
#> [4] "size_group" "mass_g" "log_mass"
#> [7] "mass_norm" "log_mass_norm" "temp_c"
#> [10] "temp_arr" "median_temp" "above_peak_temp"
#> [13] "common_name" "species" "species_ab"
#> [16] "env_temp_min" "env_temp_max" "log_y"
#> [19] "log_mass_intra" "temp_arr_intra" "log_mass_ct"
#> [22] "temp_arr_ct"
ggplot(dat, aes(temp_c, mass_g*y, color = factor(mass_g))) +
geom_point() +
facet_wrap(~species_ab, scales = "free") +
theme_classic() +
stat_smooth(se = FALSE) +
guides(color = FALSE) +
scale_color_viridis(discrete = TRUE)
ggplot(dat, aes(temp_c, y, group = factor(mass_g))) +
facet_wrap(~species_ab, scales = "free") +
theme_classic() +
stat_smooth(se = FALSE, color = "gray") +
geom_point(data = dat, aes(temp_c, y, color = above_peak_temp))
ggplot(dat, aes(temp_arr, log_y, group = factor(mass_g))) +
facet_wrap(~species_ab, scales = "free", ncol = 3) +
theme_classic() +
stat_smooth(method = "lm", se = FALSE, color = "gray") +
geom_point(data = dat, aes(temp_arr, log_y, color = above_peak_temp))
# No only using sub-optimum temperatures
ggplot(filter(dat, above_peak_temp == "N"), aes(temp_arr, log_y, group = factor(mass_g))) +
facet_wrap(~species_ab, scales = "free", ncol = 3) +
theme_classic() +
stat_smooth(method = "lm", se = FALSE, color = "gray") +
geom_point(data = filter(dat, above_peak_temp == "N"),
aes(temp_arr, log_y, color = above_peak_temp))
Fit models
- get priors
m_gauss <- brm(
log_y ~ log_mass_ct*temp_arr_ct + (1 | species_ab),
data = filter(dat, above_peak_temp == "N"),
family = gaussian(), save_pars = save_pars(all = TRUE),
control = list(adapt_delta = 0.95, max_treedepth = 12),
iter = 4000, cores = 4, chains = 4)
#> Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
#> clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Users/maxlindmark/Library/R/4.0/library/Rcpp/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/unsupported" -I"/Users/maxlindmark/Library/R/4.0/library/BH/include" -I"/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/src/" -I"/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppParallel/include/" -I"/Users/maxlindmark/Library/R/4.0/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
#> In file included from <built-in>:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Core:88:
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
#> namespace Eigen {
#> ^
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
#> namespace Eigen {
#> ^
#> ;
#> In file included from <built-in>:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#> #include <complex>
#> ^~~~~~~~~
#> 3 errors generated.
#> make: *** [foo.o] Error 1
m_student <- brm(
log_y ~ log_mass_ct*temp_arr_ct + (1 | species_ab),
data = filter(dat, above_peak_temp == "N"),
family = student(), save_pars = save_pars(all = TRUE),
control = list(adapt_delta = 0.95, max_treedepth = 12),
iter = 4000, cores = 4, chains = 4)
#> Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
#> clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Users/maxlindmark/Library/R/4.0/library/Rcpp/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/unsupported" -I"/Users/maxlindmark/Library/R/4.0/library/BH/include" -I"/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/src/" -I"/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppParallel/include/" -I"/Users/maxlindmark/Library/R/4.0/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
#> In file included from <built-in>:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Core:88:
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
#> namespace Eigen {
#> ^
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
#> namespace Eigen {
#> ^
#> ;
#> In file included from <built-in>:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#> #include <complex>
#> ^~~~~~~~~
#> 3 errors generated.
#> make: *** [foo.o] Error 1
m_skew <- brm(
log_y ~ log_mass_ct*temp_arr_ct + (1 | species_ab),
data = filter(dat, above_peak_temp == "N"),
family = skew_normal(), save_pars = save_pars(all = TRUE),
control = list(adapt_delta = 0.95, max_treedepth = 12),
iter = 4000, cores = 4, chains = 4)
#> Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
#> clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Users/maxlindmark/Library/R/4.0/library/Rcpp/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/unsupported" -I"/Users/maxlindmark/Library/R/4.0/library/BH/include" -I"/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/src/" -I"/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppParallel/include/" -I"/Users/maxlindmark/Library/R/4.0/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
#> In file included from <built-in>:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Core:88:
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
#> namespace Eigen {
#> ^
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
#> namespace Eigen {
#> ^
#> ;
#> In file included from <built-in>:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#> #include <complex>
#> ^~~~~~~~~
#> 3 errors generated.
#> make: *** [foo.o] Error 1
loo_gauss <- loo(m_gauss, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
loo_student <- loo(m_student, moment_match = TRUE)
loo_skew <- loo(m_skew, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
loo_compare(loo_gauss, loo_student, loo_skew)
#> elpd_diff se_diff
#> m_skew 0.0 0.0
#> m_student -4.8 5.8
#> m_gauss -19.5 5.6Make plots
# Summary
# https://cran.r-project.org/web/packages/sjPlot/vignettes/tab_model_estimates.html
tab_model(m_skew)| Â | log y | |
|---|---|---|
| Predictors | Estimates | CI (95%) |
| Intercept | 0.63 | 0.29 â 0.96 |
| log mass ct | -0.37 | -0.40 â -0.33 |
| temp arr ct | -0.61 | -0.71 â -0.51 |
| log_mass_ct:temp_arr_ct | 0.02 | -0.01 â 0.05 |
| N species_ab | 13 | |
| Observations | 156 | |
| Marginal R2 / Conditional R2 | 0.748 / 0.884 | |
# Plot
conditional_effects(m_skew)
ggplot(dat, aes(temp_c, temp_arr)) +
geom_point() +
stat_smooth(method = "lm") +
geom_hline(yintercept = c(39, 41), color = "tomato") +
geom_vline(xintercept = c(10, 24), color = "tomato")
pal <- brewer.pal(n = 3, name = "Set1")
summary(m_skew)
#> Family: skew_normal
#> Links: mu = identity; sigma = identity; alpha = identity
#> Formula: log_y ~ log_mass_ct * temp_arr_ct + (1 | species_ab)
#> Data: filter(dat, above_peak_temp == "N") (Number of observations: 156)
#> Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#> total post-warmup draws = 8000
#>
#> Group-Level Effects:
#> ~species_ab (Number of levels: 13)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.59 0.16 0.37 0.97 1.00 1560 2870
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 0.63 0.17 0.29 0.96 1.00 1468
#> log_mass_ct -0.37 0.02 -0.40 -0.33 1.00 5387
#> temp_arr_ct -0.61 0.05 -0.71 -0.51 1.00 4433
#> log_mass_ct:temp_arr_ct 0.02 0.02 -0.01 0.05 1.00 5485
#> Tail_ESS
#> Intercept 2034
#> log_mass_ct 4798
#> temp_arr_ct 4435
#> log_mass_ct:temp_arr_ct 5782
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 0.30 0.02 0.26 0.34 1.00 5798 5050
#> alpha -6.54 2.00 -11.02 -3.38 1.00 4767 4127
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
p1 <- data.frame(
expand.grid(log_mass_ct = seq_range(dat$log_mass_ct, n = 101),
temp_arr_ct = round(c(39, 41) - mean(dat$temp_arr), 1))) %>%
mutate(temp_arr = round(temp_arr_ct + mean(dat$temp_arr), 1)) %>%
mutate(temp_c = round(1/(temp_arr * (8.617*10^-5)) - 273.15, 0)) %>%
mutate(temp = paste(temp_arr, " (", temp_c, "°C)", sep = "")) %>%
add_predicted_draws(m_skew, re_formula = NA) %>%
ggplot(., aes(x = log_mass_ct, y = .prediction,
color = factor(temp), fill = factor(temp))) +
stat_lineribbon(.width = c(.85, .95), alpha = 1/4) +
stat_lineribbon(.width = 0, alpha = 0.7) +
scale_color_brewer(palette = "Set1", name = "") +
scale_fill_brewer(palette = "Set1", name = "Temperature") +
theme_classic(base_size = 10) +
geom_point(data = dat, aes(log_mass_ct, log_y), inherit.aes = FALSE,
size = 1.5, alpha = 0.6, shape = 21, color = "white", fill = "grey30") +
# guides(color = FALSE, fill = guide_legend(ncol = 1)) +
guides(color = FALSE) +
labs(x = "log(mass)") +
# annotate("text", 0, 4, size = 3, color = pal[2],
# label = "y=0.64 - 0.37Ălog(m) - 0.61Ă", fontface = "italic") + # Cold
# annotate("text", 0, 3.5, size = 3, color = pal[1],
# label = "y=0.64 - 0.37Ălog(m) - 0.61Ăt", fontface = "italic") + # Warm
theme(aspect.ratio = 1,
legend.key.size = unit(0.5, "line"),
legend.text = element_text(size = 6),
legend.title = element_text(size = 7),
legend.spacing.x = unit(0.1, 'cm'),
legend.margin = margin(0, 0, 0, 0),
#legend.position = c(0.3, 0.2),
legend.position = "bottom")
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
p2 <- data.frame(
expand.grid(log_mass_ct = round(c(2, 5) - mean(dat$log_mass), 2), # summary(dat$log_mass)
temp_arr_ct = seq_range(dat$temp_arr_ct, n = 101))) %>%
mutate(log_mass = round(log_mass_ct + mean(dat$log_mass), 2)) %>%
add_predicted_draws(m_skew, re_formula = NA) %>%
ggplot(., aes(x = temp_arr_ct, y = .prediction,
color = factor(log_mass), fill = factor(log_mass))) +
stat_lineribbon(.width = c(.85, .95), alpha = 1/4) +
stat_lineribbon(.width = 0, alpha = 0.7) +
scale_color_brewer(palette = "Dark2", name = "") +
scale_fill_brewer(palette = "Dark2", name = "log(mass)") + # Check this! Not the centered scale
theme_classic(base_size = 10) +
geom_point(data = dat, aes(temp_arr_ct, log_y), inherit.aes = FALSE,
size = 1.5, alpha = 0.6, shape = 21, color = "white", fill = "grey30") +
#guides(color = FALSE, fill = guide_legend(ncol = 1)) +
guides(color = FALSE) +
labs(x = "Temperature (1/kT[K])") +
theme(aspect.ratio = 1,
legend.key.size = unit(0.5, "line"),
legend.text = element_text(size = 6),
legend.title = element_text(size = 7),
legend.spacing.x = unit(0.1, 'cm'),
legend.margin = margin(0, 0, 0, 0),
#legend.position = c(0.2, 0.2),
legend.position = "bottom")
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
p3 <- m_skew %>%
spread_draws(`b_log_mass_ct:temp_arr_ct`) %>%
ggplot(aes(x = `b_log_mass_ct:temp_arr_ct`, fill = stat(x < 0))) +
stat_halfeye(.width = c(0.85, 0.95)) +
theme_classic(base_size = 10) +
geom_vline(xintercept = 0, linetype = "dashed") +
scale_fill_manual(values = c("gray90", "grey60")) +
labs(x = "MassĂTemp. interaction") +
#guides(fill = guide_legend(ncol = 1)) +
theme(aspect.ratio = 1,
legend.key.size = unit(0.5, "line"),
legend.text = element_text(size = 6),
legend.title = element_text(size = 7),
legend.spacing.x = unit(0.1, 'cm'),
legend.margin = margin(0, 0, 0, 0),
#legend.position = c(0.8, 0.5),
#legend.position = "bottom"
legend.position = "bottom")
(p1 | p2 | p3) + plot_annotation(tag_levels = "A")
#(p1 / p2 / p3) + plot_annotation(tag_levels = "A")
ggsave("figures/Fig1.png", width = 6.5, height = 6.5, dpi = 600)
get_variables(m_skew)
#> [1] "b_Intercept"
#> [2] "b_log_mass_ct"
#> [3] "b_temp_arr_ct"
#> [4] "b_log_mass_ct:temp_arr_ct"
#> [5] "sd_species_ab__Intercept"
#> [6] "sigma"
#> [7] "alpha"
#> [8] "Intercept"
#> [9] "r_species_ab[A.minor,Intercept]"
#> [10] "r_species_ab[B.saida,Intercept]"
#> [11] "r_species_ab[C.lumpus,Intercept]"
#> [12] "r_species_ab[G.morhua,Intercept]"
#> [13] "r_species_ab[H.hippoglossus,Intercept]"
#> [14] "r_species_ab[L.calcarifer,Intercept]"
#> [15] "r_species_ab[P.fulvidraco,Intercept]"
#> [16] "r_species_ab[P.olivaceus,Intercept]"
#> [17] "r_species_ab[P.yokohamae,Intercept]"
#> [18] "r_species_ab[R.canadum,Intercept]"
#> [19] "r_species_ab[S.alpinus,Intercept]"
#> [20] "r_species_ab[S.maximus,Intercept]"
#> [21] "r_species_ab[S.salar,Intercept]"
#> [22] "lp__"
#> [23] "z_1[1,1]"
#> [24] "z_1[1,2]"
#> [25] "z_1[1,3]"
#> [26] "z_1[1,4]"
#> [27] "z_1[1,5]"
#> [28] "z_1[1,6]"
#> [29] "z_1[1,7]"
#> [30] "z_1[1,8]"
#> [31] "z_1[1,9]"
#> [32] "z_1[1,10]"
#> [33] "z_1[1,11]"
#> [34] "z_1[1,12]"
#> [35] "z_1[1,13]"
#> [36] "accept_stat__"
#> [37] "stepsize__"
#> [38] "treedepth__"
#> [39] "n_leapfrog__"
#> [40] "divergent__"
#> [41] "energy__"
p4 <- m_skew %>%
spread_draws(b_Intercept, r_species_ab[species_ab, ]) %>%
mutate(intercept_mean = b_Intercept + r_species_ab) %>%
ggplot(aes(y = species_ab, x = intercept_mean)) +
stat_halfeye() +
theme_classic(base_size = 12) +
#coord_cartesian(xlim = c(-0.2, 0.25)) +
scale_fill_manual(values = c("gray80", "tomato")) +
theme(legend.position = "bottom",
axis.text.y = element_text(face = "italic")) +
labs(y = "Species", x = "Intercept") +
NULL
p4
ggsave("figures/Fig2.png", width = 3, height = 6, dpi = 600)
# Model validation
posterior <- as.array(m_skew)
dimnames(posterior)
#> $iteration
#> [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10"
#> [11] "11" "12" "13" "14" "15" "16" "17" "18" "19" "20"
#> [21] "21" "22" "23" "24" "25" "26" "27" "28" "29" "30"
#> [31] "31" "32" "33" "34" "35" "36" "37" "38" "39" "40"
#> [41] "41" "42" "43" "44" "45" "46" "47" "48" "49" "50"
#> [51] "51" "52" "53" "54" "55" "56" "57" "58" "59" "60"
#> [61] "61" "62" "63" "64" "65" "66" "67" "68" "69" "70"
#> [71] "71" "72" "73" "74" "75" "76" "77" "78" "79" "80"
#> [81] "81" "82" "83" "84" "85" "86" "87" "88" "89" "90"
#> [91] "91" "92" "93" "94" "95" "96" "97" "98" "99" "100"
#> [101] "101" "102" "103" "104" "105" "106" "107" "108" "109" "110"
#> [111] "111" "112" "113" "114" "115" "116" "117" "118" "119" "120"
#> [121] "121" "122" "123" "124" "125" "126" "127" "128" "129" "130"
#> [131] "131" "132" "133" "134" "135" "136" "137" "138" "139" "140"
#> [141] "141" "142" "143" "144" "145" "146" "147" "148" "149" "150"
#> [151] "151" "152" "153" "154" "155" "156" "157" "158" "159" "160"
#> [161] "161" "162" "163" "164" "165" "166" "167" "168" "169" "170"
#> [171] "171" "172" "173" "174" "175" "176" "177" "178" "179" "180"
#> [181] "181" "182" "183" "184" "185" "186" "187" "188" "189" "190"
#> [191] "191" "192" "193" "194" "195" "196" "197" "198" "199" "200"
#> [201] "201" "202" "203" "204" "205" "206" "207" "208" "209" "210"
#> [211] "211" "212" "213" "214" "215" "216" "217" "218" "219" "220"
#> [221] "221" "222" "223" "224" "225" "226" "227" "228" "229" "230"
#> [231] "231" "232" "233" "234" "235" "236" "237" "238" "239" "240"
#> [241] "241" "242" "243" "244" "245" "246" "247" "248" "249" "250"
#> [251] "251" "252" "253" "254" "255" "256" "257" "258" "259" "260"
#> [261] "261" "262" "263" "264" "265" "266" "267" "268" "269" "270"
#> [271] "271" "272" "273" "274" "275" "276" "277" "278" "279" "280"
#> [281] "281" "282" "283" "284" "285" "286" "287" "288" "289" "290"
#> [291] "291" "292" "293" "294" "295" "296" "297" "298" "299" "300"
#> [301] "301" "302" "303" "304" "305" "306" "307" "308" "309" "310"
#> [311] "311" "312" "313" "314" "315" "316" "317" "318" "319" "320"
#> [321] "321" "322" "323" "324" "325" "326" "327" "328" "329" "330"
#> [331] "331" "332" "333" "334" "335" "336" "337" "338" "339" "340"
#> [341] "341" "342" "343" "344" "345" "346" "347" "348" "349" "350"
#> [351] "351" "352" "353" "354" "355" "356" "357" "358" "359" "360"
#> [361] "361" "362" "363" "364" "365" "366" "367" "368" "369" "370"
#> [371] "371" "372" "373" "374" "375" "376" "377" "378" "379" "380"
#> [381] "381" "382" "383" "384" "385" "386" "387" "388" "389" "390"
#> [391] "391" "392" "393" "394" "395" "396" "397" "398" "399" "400"
#> [401] "401" "402" "403" "404" "405" "406" "407" "408" "409" "410"
#> [411] "411" "412" "413" "414" "415" "416" "417" "418" "419" "420"
#> [421] "421" "422" "423" "424" "425" "426" "427" "428" "429" "430"
#> [431] "431" "432" "433" "434" "435" "436" "437" "438" "439" "440"
#> [441] "441" "442" "443" "444" "445" "446" "447" "448" "449" "450"
#> [451] "451" "452" "453" "454" "455" "456" "457" "458" "459" "460"
#> [461] "461" "462" "463" "464" "465" "466" "467" "468" "469" "470"
#> [471] "471" "472" "473" "474" "475" "476" "477" "478" "479" "480"
#> [481] "481" "482" "483" "484" "485" "486" "487" "488" "489" "490"
#> [491] "491" "492" "493" "494" "495" "496" "497" "498" "499" "500"
#> [501] "501" "502" "503" "504" "505" "506" "507" "508" "509" "510"
#> [511] "511" "512" "513" "514" "515" "516" "517" "518" "519" "520"
#> [521] "521" "522" "523" "524" "525" "526" "527" "528" "529" "530"
#> [531] "531" "532" "533" "534" "535" "536" "537" "538" "539" "540"
#> [541] "541" "542" "543" "544" "545" "546" "547" "548" "549" "550"
#> [551] "551" "552" "553" "554" "555" "556" "557" "558" "559" "560"
#> [561] "561" "562" "563" "564" "565" "566" "567" "568" "569" "570"
#> [571] "571" "572" "573" "574" "575" "576" "577" "578" "579" "580"
#> [581] "581" "582" "583" "584" "585" "586" "587" "588" "589" "590"
#> [591] "591" "592" "593" "594" "595" "596" "597" "598" "599" "600"
#> [601] "601" "602" "603" "604" "605" "606" "607" "608" "609" "610"
#> [611] "611" "612" "613" "614" "615" "616" "617" "618" "619" "620"
#> [621] "621" "622" "623" "624" "625" "626" "627" "628" "629" "630"
#> [631] "631" "632" "633" "634" "635" "636" "637" "638" "639" "640"
#> [641] "641" "642" "643" "644" "645" "646" "647" "648" "649" "650"
#> [651] "651" "652" "653" "654" "655" "656" "657" "658" "659" "660"
#> [661] "661" "662" "663" "664" "665" "666" "667" "668" "669" "670"
#> [671] "671" "672" "673" "674" "675" "676" "677" "678" "679" "680"
#> [681] "681" "682" "683" "684" "685" "686" "687" "688" "689" "690"
#> [691] "691" "692" "693" "694" "695" "696" "697" "698" "699" "700"
#> [701] "701" "702" "703" "704" "705" "706" "707" "708" "709" "710"
#> [711] "711" "712" "713" "714" "715" "716" "717" "718" "719" "720"
#> [721] "721" "722" "723" "724" "725" "726" "727" "728" "729" "730"
#> [731] "731" "732" "733" "734" "735" "736" "737" "738" "739" "740"
#> [741] "741" "742" "743" "744" "745" "746" "747" "748" "749" "750"
#> [751] "751" "752" "753" "754" "755" "756" "757" "758" "759" "760"
#> [761] "761" "762" "763" "764" "765" "766" "767" "768" "769" "770"
#> [771] "771" "772" "773" "774" "775" "776" "777" "778" "779" "780"
#> [781] "781" "782" "783" "784" "785" "786" "787" "788" "789" "790"
#> [791] "791" "792" "793" "794" "795" "796" "797" "798" "799" "800"
#> [801] "801" "802" "803" "804" "805" "806" "807" "808" "809" "810"
#> [811] "811" "812" "813" "814" "815" "816" "817" "818" "819" "820"
#> [821] "821" "822" "823" "824" "825" "826" "827" "828" "829" "830"
#> [831] "831" "832" "833" "834" "835" "836" "837" "838" "839" "840"
#> [841] "841" "842" "843" "844" "845" "846" "847" "848" "849" "850"
#> [851] "851" "852" "853" "854" "855" "856" "857" "858" "859" "860"
#> [861] "861" "862" "863" "864" "865" "866" "867" "868" "869" "870"
#> [871] "871" "872" "873" "874" "875" "876" "877" "878" "879" "880"
#> [881] "881" "882" "883" "884" "885" "886" "887" "888" "889" "890"
#> [891] "891" "892" "893" "894" "895" "896" "897" "898" "899" "900"
#> [901] "901" "902" "903" "904" "905" "906" "907" "908" "909" "910"
#> [911] "911" "912" "913" "914" "915" "916" "917" "918" "919" "920"
#> [921] "921" "922" "923" "924" "925" "926" "927" "928" "929" "930"
#> [931] "931" "932" "933" "934" "935" "936" "937" "938" "939" "940"
#> [941] "941" "942" "943" "944" "945" "946" "947" "948" "949" "950"
#> [951] "951" "952" "953" "954" "955" "956" "957" "958" "959" "960"
#> [961] "961" "962" "963" "964" "965" "966" "967" "968" "969" "970"
#> [971] "971" "972" "973" "974" "975" "976" "977" "978" "979" "980"
#> [981] "981" "982" "983" "984" "985" "986" "987" "988" "989" "990"
#> [991] "991" "992" "993" "994" "995" "996" "997" "998" "999" "1000"
#> [1001] "1001" "1002" "1003" "1004" "1005" "1006" "1007" "1008" "1009" "1010"
#> [1011] "1011" "1012" "1013" "1014" "1015" "1016" "1017" "1018" "1019" "1020"
#> [1021] "1021" "1022" "1023" "1024" "1025" "1026" "1027" "1028" "1029" "1030"
#> [1031] "1031" "1032" "1033" "1034" "1035" "1036" "1037" "1038" "1039" "1040"
#> [1041] "1041" "1042" "1043" "1044" "1045" "1046" "1047" "1048" "1049" "1050"
#> [1051] "1051" "1052" "1053" "1054" "1055" "1056" "1057" "1058" "1059" "1060"
#> [1061] "1061" "1062" "1063" "1064" "1065" "1066" "1067" "1068" "1069" "1070"
#> [1071] "1071" "1072" "1073" "1074" "1075" "1076" "1077" "1078" "1079" "1080"
#> [1081] "1081" "1082" "1083" "1084" "1085" "1086" "1087" "1088" "1089" "1090"
#> [1091] "1091" "1092" "1093" "1094" "1095" "1096" "1097" "1098" "1099" "1100"
#> [1101] "1101" "1102" "1103" "1104" "1105" "1106" "1107" "1108" "1109" "1110"
#> [1111] "1111" "1112" "1113" "1114" "1115" "1116" "1117" "1118" "1119" "1120"
#> [1121] "1121" "1122" "1123" "1124" "1125" "1126" "1127" "1128" "1129" "1130"
#> [1131] "1131" "1132" "1133" "1134" "1135" "1136" "1137" "1138" "1139" "1140"
#> [1141] "1141" "1142" "1143" "1144" "1145" "1146" "1147" "1148" "1149" "1150"
#> [1151] "1151" "1152" "1153" "1154" "1155" "1156" "1157" "1158" "1159" "1160"
#> [1161] "1161" "1162" "1163" "1164" "1165" "1166" "1167" "1168" "1169" "1170"
#> [1171] "1171" "1172" "1173" "1174" "1175" "1176" "1177" "1178" "1179" "1180"
#> [1181] "1181" "1182" "1183" "1184" "1185" "1186" "1187" "1188" "1189" "1190"
#> [1191] "1191" "1192" "1193" "1194" "1195" "1196" "1197" "1198" "1199" "1200"
#> [1201] "1201" "1202" "1203" "1204" "1205" "1206" "1207" "1208" "1209" "1210"
#> [1211] "1211" "1212" "1213" "1214" "1215" "1216" "1217" "1218" "1219" "1220"
#> [1221] "1221" "1222" "1223" "1224" "1225" "1226" "1227" "1228" "1229" "1230"
#> [1231] "1231" "1232" "1233" "1234" "1235" "1236" "1237" "1238" "1239" "1240"
#> [1241] "1241" "1242" "1243" "1244" "1245" "1246" "1247" "1248" "1249" "1250"
#> [1251] "1251" "1252" "1253" "1254" "1255" "1256" "1257" "1258" "1259" "1260"
#> [1261] "1261" "1262" "1263" "1264" "1265" "1266" "1267" "1268" "1269" "1270"
#> [1271] "1271" "1272" "1273" "1274" "1275" "1276" "1277" "1278" "1279" "1280"
#> [1281] "1281" "1282" "1283" "1284" "1285" "1286" "1287" "1288" "1289" "1290"
#> [1291] "1291" "1292" "1293" "1294" "1295" "1296" "1297" "1298" "1299" "1300"
#> [1301] "1301" "1302" "1303" "1304" "1305" "1306" "1307" "1308" "1309" "1310"
#> [1311] "1311" "1312" "1313" "1314" "1315" "1316" "1317" "1318" "1319" "1320"
#> [1321] "1321" "1322" "1323" "1324" "1325" "1326" "1327" "1328" "1329" "1330"
#> [1331] "1331" "1332" "1333" "1334" "1335" "1336" "1337" "1338" "1339" "1340"
#> [1341] "1341" "1342" "1343" "1344" "1345" "1346" "1347" "1348" "1349" "1350"
#> [1351] "1351" "1352" "1353" "1354" "1355" "1356" "1357" "1358" "1359" "1360"
#> [1361] "1361" "1362" "1363" "1364" "1365" "1366" "1367" "1368" "1369" "1370"
#> [1371] "1371" "1372" "1373" "1374" "1375" "1376" "1377" "1378" "1379" "1380"
#> [1381] "1381" "1382" "1383" "1384" "1385" "1386" "1387" "1388" "1389" "1390"
#> [1391] "1391" "1392" "1393" "1394" "1395" "1396" "1397" "1398" "1399" "1400"
#> [1401] "1401" "1402" "1403" "1404" "1405" "1406" "1407" "1408" "1409" "1410"
#> [1411] "1411" "1412" "1413" "1414" "1415" "1416" "1417" "1418" "1419" "1420"
#> [1421] "1421" "1422" "1423" "1424" "1425" "1426" "1427" "1428" "1429" "1430"
#> [1431] "1431" "1432" "1433" "1434" "1435" "1436" "1437" "1438" "1439" "1440"
#> [1441] "1441" "1442" "1443" "1444" "1445" "1446" "1447" "1448" "1449" "1450"
#> [1451] "1451" "1452" "1453" "1454" "1455" "1456" "1457" "1458" "1459" "1460"
#> [1461] "1461" "1462" "1463" "1464" "1465" "1466" "1467" "1468" "1469" "1470"
#> [1471] "1471" "1472" "1473" "1474" "1475" "1476" "1477" "1478" "1479" "1480"
#> [1481] "1481" "1482" "1483" "1484" "1485" "1486" "1487" "1488" "1489" "1490"
#> [1491] "1491" "1492" "1493" "1494" "1495" "1496" "1497" "1498" "1499" "1500"
#> [1501] "1501" "1502" "1503" "1504" "1505" "1506" "1507" "1508" "1509" "1510"
#> [1511] "1511" "1512" "1513" "1514" "1515" "1516" "1517" "1518" "1519" "1520"
#> [1521] "1521" "1522" "1523" "1524" "1525" "1526" "1527" "1528" "1529" "1530"
#> [1531] "1531" "1532" "1533" "1534" "1535" "1536" "1537" "1538" "1539" "1540"
#> [1541] "1541" "1542" "1543" "1544" "1545" "1546" "1547" "1548" "1549" "1550"
#> [1551] "1551" "1552" "1553" "1554" "1555" "1556" "1557" "1558" "1559" "1560"
#> [1561] "1561" "1562" "1563" "1564" "1565" "1566" "1567" "1568" "1569" "1570"
#> [1571] "1571" "1572" "1573" "1574" "1575" "1576" "1577" "1578" "1579" "1580"
#> [1581] "1581" "1582" "1583" "1584" "1585" "1586" "1587" "1588" "1589" "1590"
#> [1591] "1591" "1592" "1593" "1594" "1595" "1596" "1597" "1598" "1599" "1600"
#> [1601] "1601" "1602" "1603" "1604" "1605" "1606" "1607" "1608" "1609" "1610"
#> [1611] "1611" "1612" "1613" "1614" "1615" "1616" "1617" "1618" "1619" "1620"
#> [1621] "1621" "1622" "1623" "1624" "1625" "1626" "1627" "1628" "1629" "1630"
#> [1631] "1631" "1632" "1633" "1634" "1635" "1636" "1637" "1638" "1639" "1640"
#> [1641] "1641" "1642" "1643" "1644" "1645" "1646" "1647" "1648" "1649" "1650"
#> [1651] "1651" "1652" "1653" "1654" "1655" "1656" "1657" "1658" "1659" "1660"
#> [1661] "1661" "1662" "1663" "1664" "1665" "1666" "1667" "1668" "1669" "1670"
#> [1671] "1671" "1672" "1673" "1674" "1675" "1676" "1677" "1678" "1679" "1680"
#> [1681] "1681" "1682" "1683" "1684" "1685" "1686" "1687" "1688" "1689" "1690"
#> [1691] "1691" "1692" "1693" "1694" "1695" "1696" "1697" "1698" "1699" "1700"
#> [1701] "1701" "1702" "1703" "1704" "1705" "1706" "1707" "1708" "1709" "1710"
#> [1711] "1711" "1712" "1713" "1714" "1715" "1716" "1717" "1718" "1719" "1720"
#> [1721] "1721" "1722" "1723" "1724" "1725" "1726" "1727" "1728" "1729" "1730"
#> [1731] "1731" "1732" "1733" "1734" "1735" "1736" "1737" "1738" "1739" "1740"
#> [1741] "1741" "1742" "1743" "1744" "1745" "1746" "1747" "1748" "1749" "1750"
#> [1751] "1751" "1752" "1753" "1754" "1755" "1756" "1757" "1758" "1759" "1760"
#> [1761] "1761" "1762" "1763" "1764" "1765" "1766" "1767" "1768" "1769" "1770"
#> [1771] "1771" "1772" "1773" "1774" "1775" "1776" "1777" "1778" "1779" "1780"
#> [1781] "1781" "1782" "1783" "1784" "1785" "1786" "1787" "1788" "1789" "1790"
#> [1791] "1791" "1792" "1793" "1794" "1795" "1796" "1797" "1798" "1799" "1800"
#> [1801] "1801" "1802" "1803" "1804" "1805" "1806" "1807" "1808" "1809" "1810"
#> [1811] "1811" "1812" "1813" "1814" "1815" "1816" "1817" "1818" "1819" "1820"
#> [1821] "1821" "1822" "1823" "1824" "1825" "1826" "1827" "1828" "1829" "1830"
#> [1831] "1831" "1832" "1833" "1834" "1835" "1836" "1837" "1838" "1839" "1840"
#> [1841] "1841" "1842" "1843" "1844" "1845" "1846" "1847" "1848" "1849" "1850"
#> [1851] "1851" "1852" "1853" "1854" "1855" "1856" "1857" "1858" "1859" "1860"
#> [1861] "1861" "1862" "1863" "1864" "1865" "1866" "1867" "1868" "1869" "1870"
#> [1871] "1871" "1872" "1873" "1874" "1875" "1876" "1877" "1878" "1879" "1880"
#> [1881] "1881" "1882" "1883" "1884" "1885" "1886" "1887" "1888" "1889" "1890"
#> [1891] "1891" "1892" "1893" "1894" "1895" "1896" "1897" "1898" "1899" "1900"
#> [1901] "1901" "1902" "1903" "1904" "1905" "1906" "1907" "1908" "1909" "1910"
#> [1911] "1911" "1912" "1913" "1914" "1915" "1916" "1917" "1918" "1919" "1920"
#> [1921] "1921" "1922" "1923" "1924" "1925" "1926" "1927" "1928" "1929" "1930"
#> [1931] "1931" "1932" "1933" "1934" "1935" "1936" "1937" "1938" "1939" "1940"
#> [1941] "1941" "1942" "1943" "1944" "1945" "1946" "1947" "1948" "1949" "1950"
#> [1951] "1951" "1952" "1953" "1954" "1955" "1956" "1957" "1958" "1959" "1960"
#> [1961] "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969" "1970"
#> [1971] "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979" "1980"
#> [1981] "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989" "1990"
#> [1991] "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999" "2000"
#>
#> $chain
#> [1] "1" "2" "3" "4"
#>
#> $variable
#> [1] "b_Intercept"
#> [2] "b_log_mass_ct"
#> [3] "b_temp_arr_ct"
#> [4] "b_log_mass_ct:temp_arr_ct"
#> [5] "sd_species_ab__Intercept"
#> [6] "sigma"
#> [7] "alpha"
#> [8] "Intercept"
#> [9] "r_species_ab[A.minor,Intercept]"
#> [10] "r_species_ab[B.saida,Intercept]"
#> [11] "r_species_ab[C.lumpus,Intercept]"
#> [12] "r_species_ab[G.morhua,Intercept]"
#> [13] "r_species_ab[H.hippoglossus,Intercept]"
#> [14] "r_species_ab[L.calcarifer,Intercept]"
#> [15] "r_species_ab[P.fulvidraco,Intercept]"
#> [16] "r_species_ab[P.olivaceus,Intercept]"
#> [17] "r_species_ab[P.yokohamae,Intercept]"
#> [18] "r_species_ab[R.canadum,Intercept]"
#> [19] "r_species_ab[S.alpinus,Intercept]"
#> [20] "r_species_ab[S.maximus,Intercept]"
#> [21] "r_species_ab[S.salar,Intercept]"
#> [22] "lp__"
#> [23] "z_1[1,1]"
#> [24] "z_1[1,2]"
#> [25] "z_1[1,3]"
#> [26] "z_1[1,4]"
#> [27] "z_1[1,5]"
#> [28] "z_1[1,6]"
#> [29] "z_1[1,7]"
#> [30] "z_1[1,8]"
#> [31] "z_1[1,9]"
#> [32] "z_1[1,10]"
#> [33] "z_1[1,11]"
#> [34] "z_1[1,12]"
#> [35] "z_1[1,13]"
mcmc_trace(
posterior,
pars = c("b_Intercept", "b_log_mass_ct",
"b_temp_arr_ct", "b_log_mass_ct:temp_arr_ct", "sd_species_ab__Intercept",
"sd_species_ab__Intercept", "sigma", "alpha"),
facet_args = list(ncol = 2, strip.position = "left")) +
geom_line(alpha = .2) +
theme(text = element_text(size = 12),
strip.text = element_text(size = 3),
legend.position = c(0.6, 0.05),
legend.direction = "horizontal") +
scale_color_viridis(discrete = TRUE, direction = -1) +
theme_classic(base_size = 12) +
theme(axis.text.x = element_text(angle = 90),
#legend.position = "bottom",
legend.position = c(0.75, 0.1),
strip.text = element_text(size = 5),
legend.direction = "horizontal") +
NULL
ggsave("figures/FigS1.png", width = 6.5, height = 6.5, dpi = 600)
# Posterior predictive
pp_check(m_skew) +
theme(text = element_text(size = 12),
legend.position = c(0.15, 0.95),
legend.background = element_rect(fill = NA)) +
scale_color_brewer(palette = "Dark2") +
labs(color = "") +
theme_classic(base_size = 14) +
theme(legend.position = c(0.2, 0.9)) +
NULL
ggsave("figures/FigS2.png", width = 4, height = 4, dpi = 600)